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ABSTRACT 



15 Recent studies indicate that altimetric observations of the ocean's mesoscale eddy field reflect 

16 the combined influence of surface buoyancy and interior potential vorticity anomalies. The 

17 former have a surface-trapped structure, while the latter have a more grave form. To assess 
is the relative importance of each contribution to the signal, it is useful to project the observed 

19 field onto a set of modes that separates their influence in a natural way. However, the 

20 surface-trapped dynamics are not well-represented by standard baroclinic modes; moreover, 

21 they are dependent on horizontal scale. 

22 Here we derive a modal decomposition that results from the simultaneous diagonalization 

23 of the energy and a generalization of potential enstrophy that includes contributions from 

24 the surface buoyancy fields. This approach yields a family of orthonomal bases that depend 

25 on two parameters: the standard baroclinic modes are recovered in a limiting case, while 

26 other choices provide modes that represent surface and interior dynamics in an efficient way. 

27 For constant stratification, these modes consist of symmetric and antisymmetric expo- 

28 nential modes that capture the surface dynamics, and a series of oscillating modes that 

29 represent the interior dynamics. Motivated by the ocean, where shears are concentrated 

30 near the upper surface, we also consider the special case of a quiescent lower surface. In this 

31 case, the interior modes are independent of wavenumber, and there is a single exponential 

32 surface mode that replaces the barotropic mode. We demonstrate the use and effectiveness 

33 of these modes by projecting the energy in a set of simulations of baroclinic turbulence. 



34 1. Introduction 



Because direct observations of the ocean's interior are sparse, satellite altimetry plays a 
crucial role in determining its time- dependent, three-dimensional velocity structure. This 
indirect measurement process assumes that sea surface height variations are dominated by 
currents with low-mode vertical structure, a result of the stiffening action of rotation and 
ensuing barotropization. Observations provide some support for this assumption, at least 
on lateral scales of order the first internal deformation scale and above. For example, using 



currentmeter records in conjunction with satellite obervations, Wunsch (1997) argues that 



the bulk of the ocean's eddy kinetic energy resides in the barotropic and first baroclinic 
modes. In addition, a number of studies show a strong correlation between the lateral size 



of eddies and the first internal deformation scale (e.g. Stammer 1997; Chelton et al. 2011). 

However, recent theoretical developments, supported by simulation and improved analysis 
of satellite altimetry, suggest that surface signals are not well-correlated with low-mode 



vertical structure, especially for submesoscale motions. In particular, Lapeyre and Klein 



48 (2006) argue that surface buoyancy and upper-ocean potential vorticity are anti-correlated 



for eddying flow, and that the three-dimensional velocity field may be obtained, assuming 
quasigeostrophy, from knowledge of the surface buoyancy field alone. The dynamics at the 
upper surface in this view are closely related to the surface quasigeostrophic (SQG) model 



52 (Blumen 1982 Held et al. 1995), and imply a vertical structure with a surface-trapped 



53 component that is not well represented by standard baroclinic modes. This view is supported 



54 by results from simulations (LaCasce and Mahadevan 2006 Klein et al. 2008), as well as 



recent analyses of satellite altimetry (e.g. Isern-Fontanet et al. 2006 Le Traon et al. 2008) 



Finally, in an atmospheric context, Tulloch and Smith (2009) have shown that lateral surface 



buoyancy gradients may interact with interior mean potential vorticity gradients to excite 
baroclinically unstable modes that generate SQG-like dynamics near the upper surface. In 
simulations, the resulting kinetic energy spectrum near the surface exhibits a steep —3 slope 
just below the deformation scale, and a flatter —5/3 slope at smaller scales — translated to 



1 



ei the oceanic context, this implies an energetic submesoscale dominated by the surface mode. 

62 One of the most widely used tools in oceanography is the projection of the vertical struc- 

63 ture of observed or simulated currents on simple bases of functions. The above observations 

64 and modeling results lead one to seek projection bases that faithfully represent both the low- 

65 mode interior structure and the surface dynamics. The standard basis of baroclinic modes, 
ee consisting of the eigenfunctions (f>(z) of the operator d z [f 2 /N 2 (z) d z (j)}, with homogenous 
67 boundary conditions d z (p\ z= o = d z (p\ z= -H = 0, fails in this respect. By construction, it is a 
es complete basis in which to expand the streamfunction ip of flows provided they satisfy the 

69 same homogeneous boundary conditions, which imply zero surface and bottom buoyancy. 

70 But for realistic flows with non-zero surface buoyancy b = fd z i/j\ z= o, expansion in baroclinic 

71 modes leads to a non-uniform convergence near z = 0, and a very large set of modes is 

72 required to capture the near-surface behaviour. 



73 As noted by Lapeyre and Klein (2006), in quasigeostrophic theory, the dynamical contri- 

74 bution of the surface buoyancy can be separated from that of the interior potential vorticity: 

75 taking advantage of the linearity of the inversion of the quasigeostrophic potential vorticity 

76 (PV) 

q = V 2 iJ + d z ({- d z A (1) 



^N 2 

77 the streamfunction may be decomposed into interior and surface parts, ip = ip mt + ^ surf 

78 (assuming zero buoyancy at the bottom), where ip mt satisfies Q with boundary condition 

79 d z ij int \ z=0 = while ^ surf satisfies the zero-PV condition VV urf + d z (f 2 /N 2 <9^ surf ) = 
so with d z ip smi \ z= o = b/f. The vertical structure of the interior contribution can be expanded 
si in the standard baroclinic modes. By contrast, the surface contribution — the only one 

82 retained in SQG theory — has a vertical structure determined by the zero PV condition 

83 which couples horizontal and vertical dependence, reducing to exp(nNz/ /), where k is the 

84 horizontal wavenumber, in the case of constant N and for z 3> H. 

85 It is intuitively clear that an effective projection basis should somehow combine modes 

86 similar to the baroclinic modes with modes that, like the exponential modes of SQG theory, 
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87 capture the dynamical contribution of the surface buoyancy. A systematic method to obtain 



ss such a basis has remained elusive, however. Tulloch and Smith (2009) proposed a heuristic 
89 model based on a barotropic and first baroclinic mode, appended by exponential modes for 



90 each surface. Similarly, Lapeyre (2009) attempted to represent the full dynamics of the 

91 upper ocean with a truncated set of standard baroclinic modes appended by an exponential 

92 surface mode. However, these hybrid modes do not diagonalize the energy, since the surface 

93 and interior modes are not orthogonal. Moreover, because the surface modes depend on 

94 wavenumber while the interior modes do not, the energetic overlap varies with horizontal 

95 scale, increasing with increasing scale. These difficulties stem from the fact that the addition 

96 of the exponential mode makes the basis functions linearly dependent in a certain sense, 

97 leading to an overcomplete frame rather than a basis. A consequence is that the modal 



98 decomposition is non-unique. Lapeyre (2009) defined a unique basis by requiring that it 

99 minimizes a certain functional, but the results remained inconclusive. An alternative basis, 

100 involving modes satisfying the Dirichlet condition ip\ z= o = together with the barotropic 



101 mode, has recently been proposed by Scott and Furnival (2012) but this too suffers from a 

102 lack of orthogonality. 

103 In this paper, we take a different approach and propose a new modal basis (or rather 

104 a family of bases) that diagonalizes the energy and effectively captures surface-intensified 

105 motion driven by buoyancy. Our approach relies on the observation that there are infinitely 
loe many possible (complete) bases onto which the flow may be projected which diagonalize the 
107 energy. As we show, a useful basis is obtained by demanding that it simultaneously diago- 
los nalizes both the energy and another quadratic invariant that generalizes potential enstrophy 
109 to include the variances of the surface and bottom buoyancy fields. The relative weight of 
no the potential enstrophy and buoyancy variances in this invariant provide two parameters 
in that determine the basis uniquely. 

n2 The eigenvalue problem that arises is similar to the standard vertical mode problem, but 

in retains a dependence on horizontal wavenumber, and the eigenvalue appears in both the 
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114 eigenvalue equation and its boundary conditions. In a limiting case, the standard baroclinic 

us modes are recovered — for constant iV and —H < z < 0, these are ipn cos(nnz/H), n = 

lie 0, 1, . . .. Another limiting case, motivated by the ocean where shears are concentrated near 

ii7 the upper surface but are weak at depth, leads to the simple basis 

ipo oc cosh [Nk(z + H)/f] , ip n oc sin [(n - l/2)nz/H)] , n = 1, 2, . . . (2) 

us which includes the exponential mode of SQG theory. 

n9 The paper is organized as follows. In section 2 we construct a generalized eigenvalue 

120 problem that defines the new basis. In section 3, we derive analytical solutions and general 

121 results for two special cases: constant N, for expository purposes, and an ocean-like case, 

122 in which the lower boundary is assumed quiescent, leading to (|2]). These modes are tested 

123 in section 4 on fields generated from a set of high-resolution quasigeostrophic simulations of 

124 baroclinic turbulence. Finally, we discuss and conclude in section 5. 



2. Surface-aware basis 

Throughout the paper, we assume a horizontally-periodic domain bounded vertically by 
rigid surfaces at z = z~ and z = z + , with total depth H = z + — z~ . The horizontal 
periodicity allows us to Fourier transform the equations in the horizontal plane, resulting in 
separable dynamics and ordinary differential equations for the vertical structure. (In more 
general domains, the Fourier series can be replaced by an expansion in eigenfunctions of 
the horizontal Laplacian, and the results obtained here should hold essentially unchanged.) 
The complex amplitudes of the quasigeostrophic potential vorticity (PV) q = qki(z), surface 
buoyancies (SBs) and streamfunction ip = tpki(z) are then related by 

P Y 

—p;ip' ] — K 2 ip = q, z < z < z + (3a) 
N z / 



./ 



2 



= &± z = z ± , (3b) 



N 2 H 
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126 where k = (k 2 + Z 2 ) 1 / 2 is the wavenumber magnitude, a prime indicates a z derivative, / 

127 is the Coriolis frequency and N = N(z) is the buoyancy frequency. We include the non- 
128 standard factor f 2 / (N 2 H) in our the definition of the SBs so that the SBs and PV have the 
129 same dimension (inverse time), and because it ultimately yields a more natural eigenvalue 
no problem. We have omitted the wavenumber subscript on q, b and ip and continue to do so 
131 onward, except where confusion may occur. 

The quasigeostrophic equation set has four quadratic invariants: energy, potential en- 
strophy, and the buoyancy variance at each surface. At each wavenumber k, these are 

E K =^ f Z+ (&m 2 + ^\ 2 ) dz 



2H J z - \N 2 

132 Summing each quantity over (k, I) gives the total invariant. 

133 We seek to define a complete basis that diagonalizes the energy. This can be done in 

134 infinitely many ways. Our strategy is based on the following principles: (i) we regard the 

135 energy as a functional, not of the streamfunction, but of the PV and of the SBs; (ii) we 

136 exploit standard results on the simultaneous diagonalization of quadratic forms. Principle 

137 (i) is grounded in the quasigeostrophic model, which makes it explicit that PV and SBs, 
us taken together, make up the set of dynamical variables. Thus, the contribution of the SBs 
139 to the dynamics is recognized; as a result, the bases we obtain naturally represent data with 
uo non-zero surface buoyancies. Regarding (ii), we recall a classical result from linear algebra: 

141 whereas there are infinitely many bases diagonalizing a quadratic form x T Ax, where A 

142 is a symmetric positive definite matrix, only one of these bases also diagonalizes another 



143 quadratic form x T Bx (e.g. Horn and Johnson 1990). This is simply found by solving the 



generalized eigenvalue problem Bx = XAx. An analogous result applies to linear operators 



145 (see, e.g. Goldstein 1980). Similarly, here we can define a unique basis by insisting that it 



146 diagonalizes another quadratic form in addition to the energy E R . A natural choice for this 



147 is a 'generalized potential enstrophy' that combines the remaining invariants into a single 
us quantity, 

P K = Z K + a + B + + a_£T (4) 

U9 where a± > are (nondimensional) undetermined weights, the choice of which will be 

150 discussed later. This approach yields a unique basis for fixed a±. 

151 To proceed, we require four objects: a vector structure that combines the SBs and 

152 interior PV, an inner product that operates on this vector, and two operators (analogous to 

153 the matrices A and B above) that give the energy and generalized potential enstrophy in 

154 terms of the inner product. These are defined as follows: 

155 Vector. We define the 'generalized potential vorticity vector]]] 

f b+ \ 



Q 



q{z) 



(5) 



156 Inner product. The specific choice of inner product is unimportant for the final results; 

157 we make what appears to be the simplest choice, namely 



(Qi,Q 2 > 



H 



qiq 2 dz + bfb% + b x b 2 



(6) 



where the overbar denotes a complex conjugate. 



159 Operators. With the definitions ^ and ([6]), it is a simple matter to find the linear operators 

160 £ and V such that 



E K = l -{Q,£Q) and P K = \{Q,VQ). 



(7) 



1 Notice that our Q bears a resemblance to the generalized potential vorticity of Bretherton ( 1966 ), which 
in our notation is written 



f 2 f 2 

ip' ) - k 2 V> - i^.i>'5(z - z + ) + jr^lp'5{z - z~). 



N 2 ^ J r N 2 ^ y N 2 
Our notation makes it plain that the PV and SBs are independent, a point that the use of Qb might obscure. 
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These are given by 



£Q 



and VQ 



a-b~ 



(8) 



-V(z) 

where the streamfunction ip is the solution of ([3]), given g and b ± . The first of these 
expressions is obtained after an integration by parts; the second is immediate. These 
two operators are positive definite and self-adjoint (see Appendix |5] for details). 



The basis we seek is now given by the eigenfunctions £ n of the generalized eigenvalue 
problem 

Vi n = fi 2 n ££ n , (9) 

where the eigenvalues y? n are positive for all n. To obtain an explicit form for Q, we 
define the components of £ n = [£+ , £ n (z), £~ ] T analogous to those of Q, and the scalar 
streamfunctions <f> n (z) such that ££ n = [(p n (z + ), — <f> n (z), — 4> n (z~)] T . In terms of these, the 
eigenvalue problem reads 



In view of (j3|, this implies that the <f) n satisfy 



/4 



P 



— K 



and 



-0n(Z 



P 



(10) 



±- 



at z = z 



± 



N 2-rnJ rn mm N 2 H a± 

This eigenvalue problem is a key result of the paper. Its eigenfunctions <f) n , which are purely 
real, give the form of the streamfunction corresponding to the basis eigenvectors £ n . The 



three components of these eigenvectors may be derived from the n using (10), although, as 
shown below, this is not necessary to project data onto the modes £ n . 

By construction, the eigenfunctions are orthogonal for the products (■,£•) and (■,?•). 
The choice of normalization for the eigenvectors £ n is inessential, but it is convenient to fix 



7 



178 the energy of each mode to be unity, that is, to take 

f 2 



N 2 



mVn + K H>mH>n 



dz &m.n. • 



(12) 



179 The expression in terms of m and <f) n is found by using ( 10 ) and ( 11 ) to eliminate £ m , £ n and 



lso the eigenvalues, then integrating by parts, which removes boundary terms. Correspondingly, 

181 

182 and 



N 2 



+ K 2 (f) m (j) n dz = (J^Sr, 



(13) 



1 

H 



dz + 



+ 



a_ 



H n 2 Smn- (14) 



The latter relation ( 14 ) has the advantage of involving only the undifferentiated streamfunc- 



184 tions, while the first relation (12) is independent of the eigenvalues and a±. 



185 The basis of eigenfunctions can be used to expand data: given Q or tp, we can write 

Q = ^ «n£n and V = ^2a n (j) n , (15) 

n n 

186 where the a n are amplitude coefficients that can be found using one of the orthogonality 



relations (12) or (13); for instance 

a n = (£,n,£Q) 



1 

H 



f_ 
N 2 



+ n 2 4> n ip dz. 



las The energy and generalized potential enstrophy are then simply 

E K =^^2\a n \ 2 and P K = - ^ fi 2 n \a n \ 2 , 

n n 

189 respectively. 



(16) 



Note that, even though the eigenvalue problem ( 11 ) is not of the standard Sturm-Liouville 
form, because of the presence of the eigenvalue /z 2 in the boundary conditions, the basis of 
eigenvectors can be shown to be complete in the sense that it provides a representation of 
arbitrary vectors Q that converges as the number of modes tends to oo. This is discussed 
further in Appendix |5j 



195 Lastly, note that our choice of orthogonality conditions implies slightly unfamiliar dimen- 

196 sions for the eigenfunctions. Because [q], [b ] ~ [T -1 ] and \p] ~ [L^ 1 ] (where T is time, L is 

197 length, and braces mean "dimensions of ") , ^ implies that [£] ~ [L~ 2 ][0]. The orthogonality 



198 condition (12) demands [0] ~ [L] and therefore [£] ~ [L 1 ]. In the next section, the problem 

199 will be analyzed in an appropriate nondimensional form. 



3. Structure of the surface-aware modes and special 



201 



cases 



202 The approach described above provides a family of bases parameterized by the values 

203 of a + and a_. In principle, different values can be chosen for different wavenumbers «; 

204 here, however, we restrict attention to choices of a± that are independent of k. To clarify 

205 some general properties of the new modes, we first recast the eigenvalue problem in non- 
206 dimensional form with the substitutions z !->■ Hz, k >->■ f/(N H)n and \i !->■ f/(NoH)/j l , 
207 where Nq is a typical value of iV; thus the wavenumber and eigenvalue are scaled by the 



208 approximate deformation length, NqH/J. The non-dimensional eigenvalue problem (11) 

209 then becomes 

A 2 + k 2 N 2 

(scj)' n ) = -X 2 n 4> n and s<j)' n = ±— n at z = 0, -1, where s = (17) 

a± l\ z {z) 

210 and we have defined an alternative eigenvalue A n such that 

li 2 n = K 2 + X 2 n . (18) 

211 Written in terms of A n , the eigenvalue equation takes the form of the standard vertical mode 

212 equation, but with more complicated boundary conditions. 



213 Analysis of the new eigenvalue problem (17) is complicated by its dependence on three 



214 independent parameters: k, a + and ct_. Moreover, for each choice of parameters, there is 

215 an infinite set of eigenvalues. Since the problem depends on the two weights a± in a nearly 

9 



216 equivalent way, we proceed first by setting the weights equal and defining a = a + = a_ (a 

217 case in which the weights differ will be considered in a later subsection). The nature of the 

218 eigenproblem is then largely determined by the size of the boundary condition coefficient 

219 fi^/ a: when fi n /a —> 0, the boundary conditions revert to the standard case 4>' n = at the 

220 top and bottom, while when fi n /a — > oo, the boundary conditions become n = at the top 

221 and bottom. However, more subtle possibilities arise as well, because unlike the standard 

222 vertical mode problem, A n may be imaginary (although //„ is always real). When A n is real, 

223 the modes are oscillatory, but when it is imaginary, the modes are evanescent — these can 

224 be interpreted either as surface modes or as extensions of the barotropic mode. 

225 This interpretation is suggested by examining the eigenvalue problem in two limiting 

226 regimes: 

227 k 2 <C a: modes with real A satisfy the simplified boundary condition (s<f)' n ) = ±\ n <ft n /a at 

228 z — 0, — 1 which further reduces to 4>' n — for a 1, corresponding to the standard 

229 baroclinic modesj^] These are complemented by a barotropic mode for which the first 

230 approximation A = can be refined to the purely imaginary A = ki/2/a. 

231 k 2 ^> a. In this case, almost all modes have n n = k 2 + X 2 n ^> a and hence satisfy the 

232 simplified boundary conditions ra = Oatz = O, — 1. There are two additional modes, 

233 however, for which fi^ = 0(a) and hence A ~ in. These solve 

u 2 

(scj)') - K 2 (j) n ~ with s(j)' n = ±—(f) n at z = 0, -1, (19) 

a± 

234 and can be recognized as surface modes, with zero interior PV. 

2 This approximation is not uniform in n but breaks down for highly oscillatory modes, with A„ = O(a), 
which satisfy (j)' = 0(a) ^ at z = 0, —1 and thus differ from the standard high-n baroclinic modes. 
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235 a. Analytical solutions for constant N 



236 In the special case of constant stratification, or s = 1, the eigenvalue problem (17) can 

237 be solved in closed form. Writing the solutions as 

<fi n = A cos(X n z) + B sm(X n z), 

238 where A and B are integration constants, and imposing the boundary conditions leads to 

239 an algebraic equation for A n , which may be either real or imaginary. For A 2 > 0, the 

240 characteristic equation (dropping the subscript n) is 

0+ + a_)A(A 2 + k 2 ) 



tan A 



(A 2 + k 2 ) 2 -a+a_A 2 ' 



(20) 



24i For A 2 < we define A = iX and obtain 



tanh A 



(a+ + ol-)X(k 2 - A 2 ) 
(k 2 - X 2 ) 2 + a + aSX 2 ' 



(21) 



Equations (20) and (21) are suitable for a graphical analysis. Fig. [I] shows that there are 



243 infinitely many solutions to (20) (top panel) and one or two solutions to (21) depending on 

244 a± (bottom panel; in both cases we set a = a + = a J). An important parameter is the ratio 



of the slopes of the right- and left-hand sides of (20) and (21) at A = 0, which in both cases 



246 IS 



a+ + a. 



k' 2 



When k < 1 there is only one solution to (21), and there is a solution of (20) with A < n/2 



248 On the other hand, if k > 1, there are two solutions to (21) (note that the maximum of the 



right-hand side of (21) is 1), and there may or may not be a solution of (20) for A < 7r/2P] 



The solution to (21 ) gives either a generalization of the barotropic mode, in the case of a 



25i single solution, or two modes that capture the vertical structure of the surface modes. Setting 



Note also that if > 4k , the denominator of the right-hand side of (20 1 goes to 0, but stays 



finite otherwise: the existence of a in the denominator determines whether there is a solution to ( 20 1 with 
A < 7r/2 in the case k~ 2 > 1. 
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252 a = a + = a_ , these solutions are plotted as functions of k in Fig. [2] there are two solutions 

253 when k > 1, but only one otherwise. The limiting solutions discussed in the previous section 



254 can be derived explicitly. In the limit k 2 = K 2 /(2a) <C 1, the single solution of (21) is given 



255 by A ~ K*J2/a, with eigenfunction oc 1, which can be interpreted as the barotropic mode. 

256 For k 2 ^> 1, the two solutions can be identified as surface intensified modes, one symmetric 

257 and the other antisymmetric about the center of the domain, explicitly given by 

0o oc cosh [k(z + |)] and 0i oc sinh [k(z + |)] , 

258 with eigenvalues fio/a = Ktanhft and fii/a = Kcoth/t. For k 3> 1, the eigenvalues are nearly 

259 identical, so that linear combinations of the eigenfunctions will also satisfy the eigenvalue 

260 problem — in particular, one can construct separate upper-surface and lower-surface modes. 



261 For real A, the right-hand side of (20) tends to zero for both large and small k, leading to 

262 eigenvalues A n = nir, n = 1,2. . . The eigenfunctions, however, differ in the two cases: for 

263 k < 1, they have the standard form n oc cos(mrz), but for k ^> 1, they are n cx sin(n7rz). 

264 The first four modes, for a = 1 and a range of n are plotted in Fig. [3j 

265 b. An oceanic special case 

Here we consider a case that is potentially the most relevant to the ocean, where shears 
near the surface may lead to surface-intensified modes, while the quiescent abyss may be 
more naturally represented by the standard boundary condition, 0' = at the bottom. The 
relevant limits for this case are a + < 1 and a_ — > oo, in which case the eigenvalue problem 
reduces to 

(s0' n )' = -A^0„, with (f) n \ z=0 = 0, 0nU=-i = 0, (22a) 

(s0' o )' - k 2 o = 0, with s0 o | z=o = — 0o, o |^=_i=O. (22b) 



266 to leading order in a + . The solutions n , n = 1, 2 . . . to (22a) describe interior modes, while 

267 O is the solution to (22b) with fil/a + = 0(1) and represents a zero PV, surface-intensified 

268 mode. 
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Note that the structure of the interior modes, like that of the standard baroclinic modes, is 
independent of k; the normalization of the mode energy that we have chosen however leads to 
K-dependent normalization factors. Since we concentrate on the leading-order approximation 
to the eigenvalue problem as a + — > 0, all the modes, including the surface-intensified one, are 
independent of a + and so are the normalisation factors (because the energy does not involve 
a+). Only the eigenvalue /x 2 , depends (linearly) on a + , although the approximation /Xq = 
can be made to conclude, in particular, that the surface-intensified mode has a generalized 
enstrophy which vanishes to leading order. 



277 Recently, Scott and Furnival (2012) proposed to use the eigenfunctions of (22a), forming 



what they term a Dirichet basis, in conjunction with the barotropic mode. While this set 
of functions, like that obtained by adding a surface mode to the standard baroclinic basis 



280 (Lapeyre 2009), does not diagonalize the energy, it is remarkable that this is achieved by 



the complete set of solutions of (22a) and (22b), that is, by the Dirichlet basis plus a surface 
mode. 



For constant N (or s — 1), the solutions to (22) may be computed explicitly; they are 
(p = A cosh [k(z + 1)] , A = 
1 



K sinh(2fi;) 



B sin 



n 



71 Z 



B 



(23a) 
(23b) 



284 in 



n 2 (n- 1/2) 2 + K 2 

with eigenvalues /xjj = a+KtanliK (corresponding to A ~ k — (a + /2) tanh k) and A n = 
l/2)7r with n = 1,2 . . .. Their dimensional form was given by ^ in the introduction. 



Again, note that the dependence on k of the coefficient for the interior modes is due to the 
normalization choice, but is irrelevant for the projection of data. 



287 4. Use of new basis for the projection of simulated data 

288 As a demonstration, we use the new basis to project the energy in three simulated tur- 

289 bulent flows, each generated by baroclinic instability of a fixed mean state in a horizontally- 
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290 periodic quasigeostrophic model. The numerical model is spectral in the horizontal, and 

291 finite-difference in the vertical — it is the same as used in, for example, Smith and Ferrari] 



292 (2009). Energy is dissipated by linear bottom drag, and enstrophy is removed by a highly 

293 scale-selective exponential cutoff filter ( Smith et al.||2002 ). In all cases, the model resolution 

294 is 512 x 512 x 100. 

295 We analyze results from three simulations. These first two are based on highly idealized 

296 flows, and will be used to demonstrate the fundamental structure of the basis, and how 

297 the partition of energy depends on both the nature of the flow, and on the choice of the 

298 nondimensional weights a±. The third simulation is based on a more realistic, ocean-like 

299 mean state, and is designed to explore the oceanic special case considered at the end of 

300 the last section. To project the simulated data onto the new basis, one must consider the 

301 generalized matrix eigenvalue problem that results from the particular vertical discretization 

302 used in the model. The details of the construction of the basis in this discretization are given 

303 explicitly Appendix B. 



304 a. Idealized 'interior' and 'surface' baroclinic instability simulations 

305 Both idealized flows have constant stratification s — 1, a ratio of domain scale to defor- 

306 mation scale equal to 4 and (3 = 0, but mean states that generate different types of baroclinic 

307 instability. The first simulation, is forced by an 'interior instability,' with a mean flow that 

308 projects onto the first (standard) baroclinic mode, U(z) = cos7rz. Flows of this type are 

309 unstable due to a sign change of the mean interior PV gradient, but have no mean SB gra- 

310 dients, since By oc U z \ z= q-i = — we refer to this simulation as BC1. The second flow 
3n is forced by an Eady mean state, with a linear mean shear U(z) = z, so the instability is 

312 driven by mean SB gradients = 1, resulting in energy generation near the two surfaces. 

313 The simulations are run to statistically steady state, and snapshots of the steady-state 

314 prognostic fields of each are used to compute horizontal (total) energy spectra. The upper 

315 panels of Fig. |4j display the horizontal spectra for the BC1 (left) and Eady (middle) simu- 
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316 lations for a few vertical levels z (the right-hand column plots will be discussed in the next 

317 subsection). It is immediately apparent that the energy in the BC1 simulation is spread 

318 rather evenly over depth; by contrast, the energy in the Eady simultion is largly concen- 

319 trated at the two surfaces. The panels in the middle row of Fig. [4] show the first few modes 

320 of the energy projected onto the standard basis, (f> n (z) oc cos(n7r,2), n = 1,2,... (the baro- 

321 clinic modes) and 0o oc 1 (the barotropic mode). Consistent with the ^-dependence of the 

322 energy in the upper panel, the energy in BC1 is largely captured by the barotropic and first 

323 baroclinic modes. By contrast, the energy in the Eady case seems to be distributed evenly 

324 across the barotropic and a large number of baroclinic modes, effectively demonstrating the 

325 failure of the standard modes to provide any insight into the energy partition in a case with 

326 large energy near the surfaces. 

327 The bottom panels of Fig. [4] display the energy spectra for the first few modes in the 

328 projection onto the new basis (BC1, left panel; Eady, middle panel). Anticipating that the 

329 BC1 simulation is best represented by the standard baroclinic basis (recovered from the 

330 generalized basis in the limit a± ^> 1), while the Eady simulation is best represented on the 

331 generalized basis in the limit a± <C 1, we chose a± = 10 6 for the former and a± = 10 -4 for 

332 the latter. As is apparent, the generalized basis with the appropriate weights more efficiently 

333 captures the surface energy in the Eady simulation much better than the standard basis. 

334 To quantify the choice of a±, we consider the projection of energy in both the BC1 

335 and Eady simulations with the generalized basis using weights ranging from a± = 1CT 3 to 

336 10 3 (always holding a = a + = a J) and ask, for what weights is the energy captured by 

337 the least number of modes? A simple diagnostic for this, the ratio of the energy contained 

338 in the first two modes to the total energy as a function of a, is shown in in Fig. [5j The 

339 results indicate that extreme values of a are best suited for the BC1 (a — > oo) and Eady 

340 (a — > 0) simulations, thus confirming our choice for Fig. |4j In the next section we examine 

341 a third simulation where the interior and surface contributions are more balanced, so that 

342 intermediate values of a± may be expected to be relevant. 
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343 b. A semi-realistic oceanic simulation 



344 The third simulation is driven by a mean state typical of the mid-latitude ocean. It uses 

345 an exponential mean stratification N 2 = N$ exp(z/h), so that s = exp(—z/h), with h = 0.2, 

346 intended to represent the pycnocline. The mean shear is U(z) = h(z + l — h) exp(z/h)+g(z) + 

347 C, where g(z) is the first standard baroclinic eigenfunction of the operator (sg')' = —\ 2 g, 

348 with g' = at z = 0, —1, so that U is surface-intensified with U'(0) = 1 and U'(—l) = 0. 

349 The constant C is set to ensure J ^ U{z) dz = 0. Both U(z) and N(z) are plotted in the top 

350 panel of Fig. [6j Note that U is baroclinically unstable due to both an internal sign change 

351 of the mean PV gradient, and to the interaction of the mean interior PV gradient Q y with 

352 the mean upper SB gradient £> + . Consistent with the assumptions of the ocean modes, the 

353 lower SB gradient B~ = 0. The ratio of the domain scale to the first baroclinic deformation 

354 radius (as determined by A -1 ) is 5. The nondimensional Coriolis gradient (5UqL~^ = 1.2, 

355 and energy is dissipated by a linear drag rL^U^ 1 = 0.4. The steady-state turbulent flow has 

356 a complicated vertical structure, as evidenced by the vertical slice of the PV shown in Fig. [7| 

357 The energy spectra for the flow are shown in the right panels of Fig. [4j just as for the BC1 

358 and Eady cases. The energy spectra by vertical level again indicates a very surface-intensified 

359 flow, but this time, the flow falls off from a —5/3 spectral slope to a more energetic interior 

360 than was the case for the Eady simulation. Projection onto the standard vertical modes 

361 (middle right panel) indicates a peak in the barotropic mode, but otherwise energy is spread 

362 evenly over a large number of baroclinic modes. Projection onto a generalized basis is shown 

363 in the bottom right panel. For this simulation with no buoyancy activity at the bottom, it 

364 is natural to use a basis with a_ — > oo. The maximum in the ratio of the energy in modes 

365 1 and 2 to total energy shown in Fig. [5] suggests that the value a = a + = 2 is appropriate. 

366 The first few modes of the corresponding basis are shown in the bottom panels of Fig [6} This 

367 is the basis chosen for Fig. [4], and indicates that the projection is very effective, with most 

368 of the energy captured by the surface and modified first baroclinic modes. An alternative 

369 basis is the 'oceanic' basis of section [b] which takes a + 1. The spectra obtained with this 
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370 basis (not shown) are essentially identical to those obtained for a + = 2. This suggests that 

371 the results are insensitive to the precise value of a + and that 'oceanic' basis may be a good 

372 default choice to analyse typical ocean data. 

373 5. Conclusion 

374 This paper presents a family of basis functions designed for the projection of three- 

375 dimensional ocean velocity data. The bases diagonalize both the quasigeostrophic energy and 

376 a generalization of the quasigeostrophic potential enstrophy that includes contributions from 

377 the buoyancy variances at the upper and lower surfaces. The family of bases is parameterized 

378 by the weights a± assigned to the surface buoyancy variances — the standard baroclinic 

379 modes are recovered in the limit a± — > oo, but the modes obtained in the opposite limit 

380 allow for efficient representation of the surface buoyancy variances. The bases should prove 

381 advantageous in a number of applications, from projection of observations to the derivation 

382 of highly truncated theoretical models. Their main drawback compared to the standard basis 

383 of baroclinic modes is the dependence of the modes on the wavenumber k which implies a lack 

384 of separation between the horizontal vertical structure in physical space. This drawback is 

385 unavoidable if some of the modes are to reflect the SQG contribution; it is minimised for the 

386 'oceanic' basis obtained for a + — » 0, a_ — > oo since all but one modes have a /t-independent 

387 structure. 

388 The limit a_ — > oo would seem a natural choice of generalized basis for typical ocean 

389 conditions takes because of the relative lack of buoyancy activity at the bottom. Regarding 

390 a + , an optimal value can in principle be chosen by inspecting the spectra for a range of 

391 values or by using a diagnostic such as that of Fig. |5j However, some simpler rules of thumb 

392 would be desirable. Intuitively, one might expect that the optimal values of a± are those 

393 that balance the contributions of the enstrophy Z K and of the surface-buoyancy variance 

394 in the generalized enstrophy P K = Z K + a + B+. Some support for this intuition is provided 
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395 by Fig. [8] which shows Z K , B K and their ratio as a function of k for the ocean simulation. 

396 The figure shows a ratio Z K /B+ that is around 5 for a broad range of k, roughly consistent 

397 with the value a + = 2 indicated by Fig. [5j There is, however, a peak around k = 4 and 

398 a substantial increase for k > 20, which suggest that better results could be obtained by 

399 allowing a + to depend on k. We have not explored this intriguing possibility here. 

400 As an alternative to the ratio Z K /B^, it would be useful to relate more directly the value 

401 of the weights a± most appropriate to project a flow on the large-scale characteristics of the 

402 flow. Since for flows driven by instabilities, Z K and are related to the large-scale PV and 

403 surface-buoyancy gradients Q y and By, it is plausible that the ratio Q y /By can be used as 

404 a guide for the choice of the weights. 

405 The advent of higher-resolution satellite observations, expected when the Surface Water 



406 Ocean Topography satellite becomes operational (Fu and Ferrari 2008), will improve our 

407 understanding of upper-ocean submesoscale dynamics only to the extent that we can connect 

408 surface observations with the three-dimensional structure of the flow below the surface. The 

409 basis derived and demonstrated here may prove a useful tool in this goal. 
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APPENDIX A 



Derivation details 

Here we prove a few relevant facts about the eigenvectors and eigenvalues of ([9]). First, 
we show that the operator £ is self-adjoint, e.g. ($ m ,££n) = (££m,£n)- Expanding the 
left-hand side and integrating by parts, we find 

1 f z+ - - -_ 

(£m, ££n) = Jj J -im<l>n + £+(f) n (z + ) - £ m <p n {z ), 

1 f Z+ ( P - Y 

= HJ ~^ n \N^^ m ) + K^rnKdz 

1 r z+ p 

\(j)' m + K 2 <p m (j) n dz, 



HJ Z - N 2 

417 since the expression on the penultimate line is clearly symmetric. The self-adjointness of V 

418 as well as the positive definiteness is obvious. 

419 To establish the completeness of the basis of the eigenvector £ n , we rewrite the eigenvalue 

420 problem in the standard form A$, n = Hn 2 ^, where A = V~ 1 £ is positive definite and self- 

421 adjoint. This operator is compact when acting on the Hilbert space of vectors Q with 

422 bounded norm (Q, Q). This is because it is essentially an integral operator with continuous 



423 kernel — the Green's function of the operator (scp 1 )' — k <f) (e.g. Debnath and Mikusihski 



section 



424 1998, section 4.8). The Hilbert-Schmidt theorem (Debnath and Mikusihski 1998, 

425 4.10) then applies to guarantee that every vector Q has a unique convergent expansion in 

426 terms of the £ n . 
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APPENDIX B 



428 Discrete eigenvalue problem and numerical computation 

429 of modes 

Here we construct the discrete version of the eigenvalue problem. Assuming a constant 
discrete coordinate Zj on J grid points, with Z\ — at the top, Zj = —H at the bottom, and 
a constant finite difference Az = Zj — Zj +1 , the mean stratification is Nq = (g/ p )Ap/ Az, 
where Ap = pj — p\ is the average background density jump between levels, pj = p(zj) 
is the background density, and p is the average density. The parameter s = Nq/N 2 is 
discretized as Sj = s(zj +1 / 2 ) = Ap/(p J+1 — pj), thus Sj is offset by a half space from pj. In 
this discretization, the SBs and PV are 

f p y i 

430 where 5 = Az/H and L D = N H/f. Nondimensionalizing /c h-> [L^, 1 ] k, ip ^ [-^|<^ _1 ] V* 

431 and (g, £» ± ) >->■ [T _1 ] (g, 6 =fc ) (for some timescale T), the discrete PV/SBs and streamfunction 

432 are related as 

Q = M>, 
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433 where 



/ 



A 



5s 1 —5si 

si -(s 1 + s 2 + 5 2 k 2 ) s 2 








Sj- 2 

.. 



(Sj_ 2 + + <T/t 2 ) Sj_i 



5s 



J-l 



(Bl) 



434 Defining the operators 



B 



1 ... 
5 ... 



. 



6 



... 1 



/ 



and F 



1 
-1 



\ 







-1 



(B2) 



one sees that B plays the part of the inner product, e.g. £ 2 ) £f B£ 2 and F accomplishes 
the awkard sign changes in the definition of the operator S. The energy in wavenumber k is 



E K 



5 



E 



5 



j-i 

3=2 



1 



For consistency with the theoretical development in section 2, we may also write the energy 
in terms of the vector Q = Aip, 

E K = ^Q*BFA- 1 Q = Iq*BSQ 

where the symmetry of F and B were used, and £ = FA -1 is defined to make the discrete 
version of the energy operator defined in ^ perfectly clear. 
Similarly, the generalized enstrophy in wavenumber k is 



P K = -Q*BVQ 
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442 where we define 



V 



\ 



a + 

1 

. . 

. . 



. . 

. . 

1 

a_ 



\ 



443 to make clear the analogy with the generalized enstrophy operator defined in (|8j). 

444 Now note that BS and BV are both symmetric (the former can be verified by checking 

445 that FBA is symmetric), so we can simultaneously diagonalize the two quadratic forms E K 

446 and P K by solving the generalized eigenvalue problem BV$,j = fi^BS^j or, in matrix form 

(BP)X = (BE )XM 2 

447 where X is the matrix with columns £j and M 2 has /i 2 along is its diagonal and zeros elsewhere. 

448 Solutions to this generalized eigenvalue problem obey the orthogonality relations 



X T B£X = I and X T BVX = M 2 , 



(B3) 



449 which are analogous to (12) and (13), respectively. 



450 In practice, it is more convenient to define a streamfunction eigenfunction such that 

451 A<p = so that the generalized eigenvalue problem can be rewritten as FVAcftj = fi^<pj, or 

452 in matrix form 

RPAcD = (DM 2 (B4) 

453 where <t> has <pj as its columns. In this case, the orthogonality relations become 



cD T FBAcD = I and <t> T "PBA 2 <t> = M 2 , 



(B5) 



454 where we've used the fact that F 2 = I. Finally, writing (B4) as <t> 1 (A X V 1 F)<t» = M 2 and 



455 using the first relation in (B5), we have the equivalent of (14) 



Q-'BV-'Q = M- 



(B6) 
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456 The expansion in the basis of eigenvectors <j> n of discrete data is readily expressed in 

457 terms of the matrix <t>. Denoting by ip the column vector of the streamfunction data (Fourier 

458 transformed in the horizontal) i[>(zj), the expansion reads 

V> = 4>a, (B7) 

459 where a = (aj, . . . , aj) T is the column vector of the mode amplitudes. These amplitudes are 

460 obtained from the data using the relation 

a = <t> T FBAi/>, 



461 which is deduced from (B5) and (B7). The total energy at a given wavenumber k, 

E K = ^*FBA^ = ^|a| 2 , 

462 where * denotes the complex (conjugate) transpose, is clearly the sum of the individual 

463 contributions \a n \ 2 /2 of each mode. Similarly, the generalized enstrophy, 

P, = \q*BVQ = i^TBAV = ^a*M 2 a, 

464 is the sum of the contributions /M^\a n \ 2 /2. 
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List of Figures 



1 Graphical solutions for eigenvalues with constant N for k = 1. The left panel 



shows the left and right hand sides of Eq. (20), and the right panel shows Eq. 



(21). 



Solutions to (21), with k scaled by \/2a, the cutoff separating cases with one 



or two solutions for imaginary A. 

The first four eigenfunctions <p n for the constant- N case, with a + = a_ = 100 
and k = 1,30,100. 

Energy spectra for the BC1 (left), Eady (middle) and Ocean (right) simula- 
tions. Top panels: spectra for selected vertical levels (see legend). Middle: 
spectra from fields projected onto standard vertical modes (modes 1, 2 and 
3-10 are shown). Bottom: spectra from fields projected onto new modes, with 
a + = a_ = 10 6 for the BC1 case, a+ = a_ = 10~ 4 for the Eady case and 
a, = 2, cn_ = 10 6 for the Ocean case. 



Ratio of the energy content of the first two modes to the total energy as a 
function of a = a + = a_ for the BC1 and Eady simulations, and as a function 
of a = a + (with a_ — > oo) for the Ocean simulation. 

Left: N 2 (z) and U(z) for the Ocean simulation. Middle: the surface mode 

<po(z) with a_ — > oo and a + < 1 (solid) and a + = 2 (dashed), for a range of 

wavenumbers k (see legend). The k — .1 lines are on top of each other. Right: 

The first three interior modes with a + C 1 and ct_ — > oo. 

Vertical slice of PV snapshot from the Ocean simulation. The flow has a 

complicated structure in the upper ocean, masking a more uniform flow at 

depth. 
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Enstrophy Z R and surface buoyancy variance as functions of wavenumber 
k for the Ocean simulation (lines with slopes -1 and -5/3 are included for 
reference). The ratio Z K /B^, also shown, can be used to guide the choice of 
the weight a + for an effective projection basis. 
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Fig. 4. Energy spectra for the BC1 (left), Eady (middle) and Ocean (right) simulations. Top 
panels: spectra for selected vertical levels (see legend). Middle: spectra from fields projected 
onto standard vertical modes (modes 1, 2 and 3-10 are shown). Bottom: spectra from fields 
projected onto new modes, with a+ = a_ = 10 6 for the BC1 case, a + = a_ = 1CT 4 for the 
Eady case and a + = 2, a_ = 10 6 for the Ocean case. 
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Fig. 5. Ratio of the energy content of the first two modes to the total energy as a function 
of a = a + = a_ for the BC1 and Eady simulations, and as a function of a = a + (with 
a_ oo) for the Ocean simulation. 
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N 2 (z) and U(z) V z > 




Fig. 6. Left: N 2 (z) and U(z) for the Ocean simulation. Middle: the surface mode (f>o(z) 
with a_ — > oo and a + C 1 (solid) and a + = 2 (dashed), for a range of wavenumbers k (see 
legend). The k — .1 lines are on top of each other. Right: The first three interior modes 
with a + < 1 and a_ — > oo. 
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Fig. 7. Vertical slice of PV snapshot from the Ocean simulation. The flow has a complicated 
structure in the upper ocean, masking a more uniform flow at depth. 
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FlG. 8. Enstrophy Z K and surface buoyancy variance as functions of wavenumber k 
for the Ocean simulation (lines with slopes -1 and -5/3 are included for reference). The 
ratio Z K /B+, also shown, can be used to guide the choice of the weight a + for an effective 
projection basis. 
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